{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Time-series ICESat1 & 2\n",
    "\n",
    "1. Plotting time-series ICESat 1\n",
    "2. Plotting time-series ICESat 1 & 2\n",
    "3. Plotting them using box and whisker plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting time-series ICESat 1\n",
    "\n",
    "### Load the ICESat 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from pathlib import Path\n",
    "import h5py\n",
    "import pandas as pd\n",
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "import cartopy.crs as ccrs\n",
    "import cartopy.io.img_tiles as cimgt\n",
    "import time\n",
    "from astropy.time import Time\n",
    "import datetime\n",
    "%matplotlib widget\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def glah14_to_df(filename,bounds):\n",
    "    ## sp_ex = [103.643, 104.667, 12.375, 13.287]\n",
    "    ## Bounds are [Longitude_West, Longitude_East, Latitude_South, Latitude_North]\n",
    "    f = h5py.File(filename, 'r')\n",
    "    lat = f['Data_40HZ']['Geolocation']['d_lat'][:]\n",
    "    lon = f['Data_40HZ']['Geolocation']['d_lon'][:]\n",
    "    elev = f['Data_40HZ']['Elevation_Surfaces']['d_elev'][:]\n",
    "    sec = f['Data_40HZ']['Elevation_Corrections']['d_satElevCorr'][:]\n",
    "    scf = f['Data_40HZ']['Quality']['sat_corr_flg'][:]\n",
    "    satndx = f['Data_40HZ']['Quality']['i_satNdx'][:]\n",
    "    dem = f['Data_40HZ']['Geophysical']['d_DEM_elv'][:]\n",
    "    date_gps = f['Data_40HZ']['DS_UTCTime_40'][1] + 630763213 \n",
    "    date_ymd = Time(date_gps,format='gps').datetime.strftime(\"%Y-%m-%d\")\n",
    "    \n",
    "    glah14_df = pd.DataFrame({'Latitude':lat,'Longitude':lon,'Elevation':elev,\n",
    "                            's_El_Corr':sec, 's_Corr_f':scf,'in_sat':satndx,\n",
    "                            'DEM':dem,'Date':date_ymd})\n",
    "    #### Subsetting\n",
    "    glah14_df_subset = glah14_df.loc[(glah14_df['Longitude']>=bounds[0]) \n",
    "                          & (glah14_df['Longitude']<=bounds[1])\n",
    "                          & (glah14_df['Latitude']>=bounds[2])\n",
    "                          & (glah14_df['Latitude']<=bounds[3])]\n",
    "    return glah14_df_subset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "##### load files\n",
    "## set the directory\n",
    "data_home = Path('/home/jovyan/ICESat_water_level/extraction/icesat/')\n",
    "## list them up and check them\n",
    "files= list(data_home.glob('*.H5'))\n",
    "### Spatial Bounds: \n",
    "tsl_sp_ex = [103.643, 104.667, 12.375, 13.287]\n",
    "### load the water mask\n",
    "tsl_wm = gpd.read_file('/home/jovyan/ICESat_water_level/extraction/shp/tsl_sample_dis.shp') ## water mask: Tonle Sap Lake"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#### convert to geoDF and clip it. \n",
    "is1_pd = pd.DataFrame()\n",
    "for num in range(1,len(files)):\n",
    "    gdf = gpd.GeoDataFrame(glah14_to_df(files[num],tsl_sp_ex))\n",
    "    gdf.set_geometry(\n",
    "        gpd.points_from_xy(gdf['Longitude'], gdf['Latitude']),\n",
    "        inplace=True, crs='EPSG:4326')\n",
    "    gdf_wm = gpd.clip(gdf, tsl_wm)\n",
    "    is1_pd = is1_pd.append(gdf_wm)\n",
    "    print(num)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### plotting one track"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "test2 = is1_lst[0]\n",
    "fig=plt.figure(figsize=(6,4))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(test2['Date'],test2['Elevation'],'.',markersize=0.25, label='all segements')\n",
    "h_leg=ax.legend()\n",
    "plt.title('Water Surface_ICESat1')\n",
    "ax.set_xlabel('Date')\n",
    "ax.set_ylabel('Water Surface, m')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Removing outliers and plot it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### outliers\n",
    "def out_rmv(df,var):\n",
    "    ### lower (25%)\n",
    "    lwq = df[var].quantile(0.25)\n",
    "    ### upper (75%)\n",
    "    hwq = df[var].quantile(0.75)\n",
    "\n",
    "    ### IQR\n",
    "    iqr = hwq - lwq\n",
    "\n",
    "    lw_out = lwq-1.5*iqr\n",
    "    hw_out = hwq+1.5*iqr\n",
    "\n",
    "    ###LSHoutlier > q0.75 + 1.5 ∗ IQR or LSHoutlier < q0.25 − 1.5 ∗ IQR\n",
    "    return df.loc[(df[var] >= lw_out) & (df[var] <= hw_out)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### I made a function for stacking the information\n",
    "def gl14_a_out(db,var2):\n",
    "    gl14_df = pd.DataFrame()\n",
    "    date_list = db['Date'].unique()\n",
    "    for prd in date_list:\n",
    "        s_db = db.loc[db['Date']==prd]\n",
    "        gl14_df = gl14_df.append(out_rmv(s_db,var2), ignore_index = True)\n",
    "    return gl14_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Average"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### I made a function for averaginf the height\n",
    "def gl14_mean(db,var2):\n",
    "    mean_df_list = pd.DataFrame()\n",
    "    date_list = db['Date'].unique()\n",
    "    for prd in date_list:\n",
    "        s_db = db.loc[db['Date']==prd]\n",
    "        av_height = statistics.mean(s_db[var2])\n",
    "        mean_df = pd.DataFrame({'Av_level':av_height,'Date':s_db['Date'].unique()})\n",
    "        mean_df_list = mean_df_list.append(mean_df, ignore_index = True)\n",
    "    return mean_df_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plottting time-series water level changes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "004ef58a08554c77aa7541a703f7a2a1",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "ename": "NameError",
     "evalue": "name 'icesat1' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-17-f7479cf45fe4>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mxaxis\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_major_locator\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mMaxNLocator\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m5\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_ylim\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m20\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m \u001b[0max\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0micesat1\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'Date'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0micesat1\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'Elevation'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'.'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcolor\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'orange'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmarkersize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.25\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlabel\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'all segements'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      7\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0micesat1\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'Date'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0micesat1\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'Elevation'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m'.'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcolor\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'green'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mmarkersize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.4\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlabel\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'Without outliers'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0;31m##ax.plot(av_h['Date'],av_h['Av_level'],'+', color = 'red',markersize=3, label='Averaged')\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'icesat1' is not defined"
     ]
    }
   ],
   "source": [
    "fig=plt.figure(figsize=(10,4))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.grid()\n",
    "ax.xaxis.set_major_locator(plt.MaxNLocator(5))\n",
    "ax.set_ylim([-20,0])\n",
    "ax.plot(icesat1['Date'],icesat1['Elevation'],'.', color = 'orange', markersize=0.25, label='all segements')\n",
    "ax.plot(icesat1['Date'],icesat1['Elevation'],'.', color = 'green',markersize=0.4, label='Without outliers')\n",
    "##ax.plot(av_h['Date'],av_h['Av_level'],'+', color = 'red',markersize=3, label='Averaged')\n",
    "##ax.plot(av_h['Date'],av_h['Av_level'],color = 'red',linewidth=0.25)\n",
    "h_leg=ax.legend()\n",
    "plt.title('Water Surface')\n",
    "ax.set_xlabel('Date')\n",
    "ax.set_ylabel('Water Surface, m')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Second way to remove outliers"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
